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Abstract. We present a detailed investigation of the probability density function 
(PDF) of order parameter fluctuations in the finite two-dimensional XY {2dXY) 
model. In the low temperature critical phase of this model, the PDF approaches 
a universal non-Gaussian limit distribution in the limit T —f 0. Our analysis 
resolves the question of temperature dependence of the PDF in this regime, for which 
conflicting results have been reported. We show analytically that a weak temperature 
dependence results from the inclusion of multiple loop graphs in a previously-derived 
graphical expansion. This is confirmed by numerical simulations on two controlled 
approximations to the 2dXY model: the Harmonic and "Harmonic XY" models. 
The Harmonic model has no Kosterlitz-Thouless-Berezinskii (KTB) transition and the 
PDF becomes progressively less skewed with increasing temperature until it closely 
approximates a Gaussian function above T « 47r. Near to that temperature we find 
some evidence of a phase transition, although our observations appear to exclude a 
thermodynamic singularity. 
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1. Introduction 

It is often assumed that global measures of many body systems will be normally 
distributed as a result of the central limit theorem (CLT). This requires that the 
system be separable into micro- or mesoscopically independent, individually negligible, 
elements [1]. When these criteria are not met, the CLT is not applicable and there 
is no reason to expect Gaussian behaviour. A prime example of this comes at 
equilibrium critical points where the infinite correlation length violates the condition of 
independence. A universal, non-Gaussian PDF of a global quantity is thus considered 
to be a signature of criticality and the characterization of the forms of these functions 
is one of the central problems in the study of critical phenomena [2]. 

The renormalization group assumes that the PDF of a global measure of a critical 
system is scale invariant and may thus be obtained from the appropriate critical 
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fixed point [3-5]. The statistics of critical fluctuations are therefore expected to be 
the same for all members of a given universality class, as evidenced by studies of 
Ising [6, 7] and Potts [8] models. Conversely there is no reason to assume that systems 
belonging to different universality classes will possess the same critical PDFs. However 
in recent years there has been much interest in the fact that experimental and numerical 
results from diverse systems, that belong to different universahty classes, in many 
cases approximate the distribution of the scalar magnetization in the two-dimensional 
XY {2dXY) model [9,10]. This model supports a continuous line of critical points 
constituting a low temperature critical 'phase' that is separated from the paramagnetic 
region by a Kosterlitz-Thouless-Berezinskii transition [11]. It was observed that the same 
limit distribution describes the magnetization fluctuations across all the fixed points on 
the critical line [12] despite the fact that they belong to different universality classes, as 
can be seen from the temperature dependence of the critical exponents [12, 13]. 

The practical and theoretical systems that show PDFs of a functional form similar to 
that of the XY model include power fiuctuations in steady state turbulence [9, 10, 14-18] 
and in liquid crystals undergoing electroconvective fiow [19-21], variations in river 
heights [22], resistance fiuctuations near electrical breakdown [23,24], numerous self- 
organised critical systems [10,25-27], dynamical fiuctuations in glassy systems [28], and 
models of equilibrium critical behaviour [29,30]. Furthermore, the distribution is related 
to the Fisher-Tippett-Gumbel distribution of extremal statistics [10] which has recently 
been shown to describe the width fluctuations of periodic, Gaussian, 1// noise [31]. It 
has been suggested that the apparent universality of the distribution may be explained 
by the hypothesis that the low temperature or spin wave limit of the 2dXY model 
captures the essence of critical behaviour for a number of universahty classes, with any 
differences that may occur being hidden by experimental error [9]. More generally, the 
low temperature limit of the 2dXY model maps on to the two-dimensional Edwards- 
Wilkinson model, a prototypical model of interface growth. Thus, these results also 
represent an important application of interface models to the study of fiuctuations in 
complex systems [32]. 

A team of authors that included one of us analytically derived an integral expression 
for the 2dXY magnetization PDF in the spin wave regime (sometimes called the BHP 
function after reference [9]) and indeed found it to be temperature-independent in 
the large- limit [13] . This derivation, which applies to the harmonic spin wave 
approximation to the XY model in the absence of vortices, provides a strong basis 
for discussions of the apparently universal form of the PDF. In contrast, Palma et 
al. [33] questioned the result of temperature independence. In their analysis of Monte 
Garlo simulations of the full 2dXY model they found evidence of a weak temperature 
dependence of the PDF throughout the low temperature phase. However, the full 2dXY 
model contains both vortex and anharmonic corrections to the harmonic spin wave 
model; although these are considered irrelevant perturbations, that is not to say that 
they do not infiuence the parameters of the PDF in finite systems. It is hard, therefore, 
to conclude with certainty from numerical studies of the 2dXY model, that the spin 
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wave contribution itself is temperature-dependent. The aim of the present work was 
to estabhsh an analytical basis for any possible temperature-dependence in the spin 
wave regime, and to test it numerically in a controlled manner by means of a series of 
approximations to the full 2dXY Hamiltonian. 

The essential properties of the models investigated are summarised in the following 
Section (Section 2) . Then in Section 3 we present a further analysis of the graphical 
expansion developed in Ref. [13]. Our numerical results are presented in Section 4. 
Firstly we simulate the PDF of the Harmonic model, which lacks both metastable vortex 
configurations and anharmonicities, so is directly comparable with the analytical results. 
We then simulate the Harmonic XY or HXY Model, in which vortices are re-introduced 
into the problem. The effect of different definitions of the order parameter is also 
investigated. Conclusions are drawn in Section 5. 



2. The 2dXY Model 



The extended critical region of the 2dXY model obviates the need for precision control of 
external constraints that one associates with locating an isolated critical point [7,8,34]. 
Furthermore the low temperature physics of the model is described precisely by a 
harmonic Hamiltonian [12, 35, 36] with the result that many properties are exactly 
calculable, without the need for renormalization or the scaling hypothesis. In the 
thermodynamic limit the magnetization, m, is zero for all finite temperatures as a result 
of the destruction of long range order by low energy spin waves - a manifestation of the 
Mermin- Wagner theorem [37]. However this limit is approached sufficiently slowly for 
m to remain physically relevant and experimentally observable [38,39]. 

The 2dXY model consists of planar spins on a two-dimensional square lattice with 
nearest neighbour interactions defined by the Hamiltonian 

H = -J J2 COs(^r - ^r')- (1) 

Here J is the (ferromagnetic) exchange interaction, 6r is the angle (relative to some 
arbitrary but fixed axis) made by the spin located at r and the sum runs over all 
nearest neighbour pairs of spins. We always deal with a periodic square lattice with 
sides L, such that N = L"^. 

Renormalization below the Kosterlitz-Thouless-Berezinskii transition temperature 
removes vortices from the problem [5, 40] and the model is described by the Harmonic 
model Hamiltonian [35] 



H = J 



(2) 



This presents a difficulty with regard to the periodicity of the spin variables implicit 
in (0). The average product of two spin variables defines a Green's function analogous 
to a propagator of a Euclidean free field. This should be able to take the full range of 
values between ±oo and so the spins in this model are necessarily non-periodic [35]. 
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Intermediate between the Harmonic and full 2dXY models is the HXy model [12], 
defined by 

H='^ Y.i^r-er'- 2n7if. (3) 

This is numerically simpler than the Villain model [36] as the parameter n is restricted 
to the values 0, ±1 to ensure that the contents of the brackets remain bounded by ±7r. 
As in the full 2dXY model, the spins are defined modulo 27r. 

The quantity we study is the instantaneous scalar magnetization, related to the 
vector order parameter by m = |m|. One definition makes use of the vector spins. 



1 

m = — A 




(4) 



We refer to this as the 'full' order parameter. A more analytically useful expression 
defines the magnetization in terms of the spin variables, 

m = — ^ cosV'r with ipr = Or — 0- (5) 

i V J. 

where 6 is the instantaneous magnetization direction. For unconfined spins this is the 
average, ^^r); ^ given configuration. For periodic spins it is convenient to define [12] 

5=1^. (6) 

2^r COS Ur 

which is the same as the mean for large A^. Defining ip enables one to work in a reference 
frame in which the Goldstone mode has been removed. 

At low temperatures angular differences between spins are likely to be small and 
it may be possible to truncate the cosine in (0). This defines the so-called 'linearized' 
order parameter, 

^ = i-^Ee (7) 



3. Evaluation of P(m) 

We are interested in the critical order parameter fiuctuations of finite systems which 
we discuss in terms of probability density functions (PDFs), P{m), calculated in the 
thermodynamic limit. Thus the order parameter must be correctly normalized to avoid 
the width of the PDF becoming either zero or infinite as A^ — oo. This has been 
extensively discussed elsewhere [9, 13] and requires that m be scaled by the standard 
deviation. We define 

n(^) = aP{m) z = (8) 

a 
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3.1. Evaluation of the Moments {m'P) 

A distribution function may be defined in terms of its moments as [41] 

-ixY 



Pirn) 



°° dx 
— ( 
-oo 27r 



p=0 



p\ 



Following from the definition of the order parameter in equation (jH)), 



^Tr (exp ( Y: ic^a^r 



\a=l 



where i is the imaginary unit and the trace is over all from 1 to and all at 
Using Gaussian integration {ipr is a Gaussian variable) we can derive [42] 

kl 



(mP) 



(m) 
2iV 



P oo 



E 

A;=0 



(9) 

(10) 
±1. 

fll) 



In this expression r = T/J, Gab is the Green's function propagator function 
rG{ra — Vb) = (i^rai^r^i ^) represents all possible pairs of a and h and the shorthand 
Hajtb = J2a J2b=ia has been introduced. The upper limits are included explicitly as a 
reminder that the sums run over p spins, not the entire lattice. 



3.2. A Graphical Representation of the Moments 

Equation (jllj) is exact. To proceed further one must develop a means of expressing the 
sums so that they may either be evaluated precisely, or approximated in a controlled 
manner. We follow [13] and consider a graphical representation of the moments which 
allows us to perform well defined partial summations. It is possible to interpret equation 
(fTT|) diagrammatically by letting Gab represent a line joining two distinct points a and 
6 on a sublattice of p points chosen from the original lattice. The /c*^ term in the 
argument of the trace is the sum of contributions from all possible graphs with k lines 
on a sublattice of p points chosen from a parent lattice of size N. 

Graphs with lines beginning and ending at the same point are disallowed (by 
virtue of the constraint a ^ b), and graphs with odd vertices (points with an odd 
number of connections) contribute zero after performing the trace. Contributions from 
disconnected graphs (those for which the lines do not form a continuous path, e.g. 
figure [T(c)| ) are equal to the product of their connected constituents. 

Relevant graphs belong to one of two categories: 'Single Loop Graphs' (SLGs) 
contain only points with no more than two lines attached; 'Multiple Loop Graphs' 
(MLGs) may have points with any (even) number of lines attached. It has previously 
been assumed [13] that MLGs do not contribute to the moments of m in the 
thermodynamic limit. We begin by reviewing the SLG only approach, to show how 
this yields a temperature independent form of n(z), before returning to consider the 
validity of neglecting MLGs. 
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Figure 1. Graphs in the Expansion of (m^): (a) Examples of allowed connected graphs 
with p — 2,3,4 and k = p. The graphs in (b) are disallowed as they contain either 
odd vertices or loops involving a single point, (c) An example of a disconnected graph 
(left hand side) which makes the same contribution as the product of its constituent 
connected parts, (d) A multiple loop graph. 



3.2.1. The Value of Single Loop Graphs Evaluation of (fTT|) requires knowledge of 
the value of each type of graph and the number of different ways each graph can 
be created. For connected SLGs it is always possible to rearrange the order of the 
Gab so that they form a simple closed loop. Transforming to Fourier space using 
G'(r) = {^/N)^^^Qe^'^'^G{(\) then gives an expression for the value of the trace in 

(HU, 



Tr 



O'aGabO'b 

ya=^b 



2P 



N 



N p 



ri=l rp=la=lbf^a 



X ^ e' 
G{k)2P 



, J2 e^'i-(r-r,)^(q^^ 

,iq,.(r,-r,)^(q^) 



(12) 
(13) 



•Ee- 

qk^o 

The factor G{k) = 2^^^{k — 1)! accounts for repetition of topologically identical 
graphs. The gu are numerical constants defined by gk = (l/iV)'^ I]q^o(l/7q)'^ with 
7q = G{q)^^ = (4 — 2cosga. — 2cosgy)~"'^. Evaluation of these constants is discussed in 
appendix B of reference [13]J 

Adjusting the factor G{k) to account for the extra combinatorial considerations 

I We note a correction to equation (B3) in reference [13]. The correct expression for is 



1 



gk = 



EE^(^2+y2) 



The numerical values that appear elsewhere in that paper are accurate so the error appears to be 
typographical. 
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imposed by disconnected graphs leads to [13] 

' V (14) 



-exp f-r) -^z^ 



The standard deviation is thus a = ^g2/2T (m), from which it is seen that hyperscahng 
is obeyed [10,13]. Substituting (fT^ into (jHl) and making the transformation x ^ x/a 
gives the BHP distribution [10,13], 

which is exphcitly independent of both system size and temperature. Equation (|15|). 
obtained using an SLG only approach, is the main result of [13]. 

3.3. The Effect of Introducing MLGs 

The simplest MLG consists of a sublattice with p points, only two of which are connected 
by > 2 lines (e.g. Figure [l(d)J . Restricting the trace to only this type of graph and 
again transforming to reciprocal space gives 



Tr 



where 



0/^ = 4 E G(qi)...G(q,)5('Eq.)- 



_ _ (17) 

qi^O,...,qfe^O 

The symmetry of G'(q) means that O2 = g^ and so, considering only the additional 
(MLG) contributions from the highly specific set of graphs discussed here, the second 
moment becomes 

(m^) = (m)^ (1 + g^T-" + le,r^ + ...+). (18) 

Therefore in the low temperature limit only the single loop graphs are significant as all 
MLGs correspond to higher powers of T, however, for high temperatures, contributions 
from MLGs may not be ruled out. 

Substituting G(q) = l/7q into (fTTj) . together with the approximation 7q q^ 
(representing the greater weight of the low frequency modes), yields, after a little algebra 

= {xl + y!){xl + yl) ^ • • • 

(19) 



(4-1 + yl-i) ((- Etl' x,)^ + (- Etl' Vi. 

The sum over {xi} is used to indicate a multiple sum over all members of the set 
{xi, X2, ■ ■ ■ ,Xk}- The terms in this sum are all positive and therefore cannot cancel each 
other. From this we can see that 0^ is explicitly non-zero: thus the contributions from 
MLGs may not be neglected in the thermodynamic limit, contrary to reference [13]. 
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(a) (b) 

Figure 2. Monte Carlo Simulations on the Harmonic Model, (a) L — 32: There is 
a clear variation with temperature from an approximately BHP form at T = 0.5 to 
a perfect Gaussian by the time T = 20.0. (b) Our simulations showed effectively no 
dependence on system size for N > 100. These results are for L = 16 and L — 32 at 
T/J = 0.5, 5.0. n and z are as defined in the text. 



4. Monte Carlo Simulations 

4.1. Simulations of the Harmonic Model 

The Harmonic model © describes tlie pliysics of tlie critical region of the 2dXY 
model [11,40]. Thus if n(z) is independent of T [13] the order parameter fluctuations 
in the Harmonic model should be BHP like at all temperatures. Single spin-flip 
Metropolis [43] Monte Carlo simulations were performed on Harmonic systems with 
L = 10, 12, 14 . . . 32, over a range of temperatures from T/J = 0.5 to as high as 
T/J = 50. Equilibration was over 10^ Monte Carlo steps per spin (MCS/s) with 10^ 
MCS/s in total. 

The results (figure |2(a)| ) show a clear variation of Tl{z) with T (for clarity only 
the data for L = 32 is shown). At low temperatures the distribution is excellently 
described by the BHP function. However as T increases the PDF becomes progressively 
less skewed, eventually becoming Gaussian. This is as expected in the light of the work 
presented in section 13.31 Higher powers of T enter the moment expansion as a result of 
the inclusion of MLGs. For r > 2, dominates the other moments at high T and the 
normalized cumulants tend to zero. We observed no variation in the form of the PDF 
as a function of system size at any temperature, as highlighted in figure |2(b)] 

The skewness, 73, provides a clear measure of the variation of the PDF with 
temperature. We have confirmed that 73 is independent of N for large systems (though 
minor corrections to this universality are observed for small A^) and a least squares 
fit to the numerical data (figure |3(a)|) reveals that the steady decrease in I73I can be 
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(a) (b) 

Figure 3. Temperature dependence of the skewness. (a) The Harmonic model with 
L = 16 over a wide range of temperatures, (b) The low temperature region for the 
Harmonic model (upper plot) with L = 16 (diamonds) and L = 32 (stars) and the 
}iXY model (lower plot) with L = 16, confirming that they obey approximately the 
same linear relationship in the critical region. The solid and dashed lines are: (a) a 
least squares quadratic fit; (b) linear regression curves; equations in the text. 

approximated (for L = 16) by (taking J = 1) 

73(r) ^ - 0.85 + 0.126r- 0.0048T2 (20) 

^ - 0.85 + 0.79r/-0.19r7^ (21) 

up to the point at which the skewness becomes zero at around rj ^ 2 (here rj = T/{2tt) 
is the spin wave exponent). The zero temperature value of 73 = —0.85 is shghtly higher 
than the theoretical BHP value of 73 = —0.89. We attribute this to the relatively small 
system used and show in figure |3(b)| that increasing the size to L = 32 gives a limiting 
value much closer to the theory with 73(T) fa —0.88 + 0.15T. We fully expect the true 
universality with respect to to be evident for larger systems. However, the skewness 
is relatively computationally expensive due to the need for averaging and as our results 
do confirm the evolution of 73 with T we leave the determination of the the precise form 
of 73 (T) from larger systems to another time. 

4-2. Temperature Dependence in the HXY Model 

Renormalization of the B.XY model below T^t removes all vortices and recovers the 
Harmonic model at some, generally higher, temperature. Using the renormalization 
procedure described by Jose et al. [40] we have calculated the mapping of temperature 
scales between these two models (table HJ showing that they are coincident below 
T < 0.9. 

Simulations of the HXF model in this range show a clear change in the skewness 
of the magnetization distribution as T is varied. Given that there are effectively no 
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Table 1. Equivalent Temperature Scales for the HXF and Harmonic Models with 
L ~ 16: The temperature of the KXY model is given as T. T^s and Kcs are the 
effective temperature and spin wave stiffness respectively - i.e. relating to the Harmonic 
model in which all vortices have been renormalized out. The RG expansion parameter 
y is a measure of the vortex density. 



vortices present in the system, it is concluded that the temperature dependence has the 
same origins as in the Harmonic model, namely the multiple loop graphs in the moment 
expansion. The plot of 73 (T) for the harmonic model (figure [3(a)| ) shows an essentially 
linear dependence at low temperatures. We confirm that a very similar form is observed 
for the HXy model (figure |3(b)| ), both L = 16 models adhering to 

73 ^ 0.827] - 0.85. (22) 



4-3. The Effect of the Order Parameter 

To enable direct comparison with the theoretical predictions all the simulations discussed 
above have used the cosine form of the order parameter, equation (0). It is interesting 
to consider the effect on the Harmonic model of substituting for this either the full order 
parameter, or the much discussed [13] linearized form ((7j). 

4-3.1. The Full Order Parameter High temperature (T > T^t) simulations of the 
2dXY and HXF models with the cosine order parameter lead to Maxwellian rather 
than Gaussian magnetization distributions. We suggest that this difference from the 
Harmonic model is not a consequence of vortices but arises from the constraints placed 
upon the spin degrees of freedom. The periodicity of the spins forces m to be positive 
and the PDF is a slice through a two-dimensional Gaussian, appearing as a Maxwell 
distribution of speeds. A similar result may be obtained for the harmonic model by 
using the full order parameter (jlj). This quantity is necessarily positive and thus has 
the same effect as constraining the spin variables to within the positive range of the 
cosine function in Monte Carlo simulations confirm this, showing BHP fluctuations 
changing to Maxwellian behaviour as T is increased (figure |3)). Despite the Harmonic 
model/full order parameter combination giving the same limiting T PDFs as the 2dXY 
and HXy models, it should be noted that the path of intermediate distributions is 
different and only for the Harmonic model is there a region of Gaussian behaviour. 
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Figure 4. Monte Carlo Results for the Harmonic Model with Full Order Parameter 
(L — 16): Again there is a clear variation of the distribution as a function of 
temperature. At low T the curve closely approximates the BHP form, but with the 
full order parameter it becomes Maxwellian, rather than Gaussian, at high T (c.f. 
figure 2(a) I. 



4-3.2. The Linearized Order Parameter The linearized order parameter ((7j) only 
describes m at temperatures low enough to exclude large angular differences between 
neighbouring spins. Despite this, Bramwell et al. demonstrated that, for the Harmonic 
model, the distribution of the linearized order parameter is precisely the BHP function 
at all temperatures [13]. We have simulated this model for a system with = 1024 
and the results confirm that the magnetization distribution is consistently of the BHP 
form for all our simulations, even up to T/J = 50. This result may now be understood 
as a manifestation of the low temperature approximation implicit in the derivation of 
the BHP distribution (|T5|) [13]. The neglect of MLGs leads to a PDF which perfectly 
describes critical fluctuations in the 2dXY model, but strictly only in the limit T ^ 0. 
The linearized OP is only a valid description of the magnetization in the same limit 
and hence has the same temperature independent fluctuations as the full OP at T = 0. 
We emphasize that we have not demonstrated an analytical equivalence between the 
neglect of MLGs and the neglect of anharmonic terms in the cosine OP. Despite this 
it seems likely that this heuristic justification of the rigorous equivalence between the 
PDFs obtained in each case is correct. 

5. Is There A Phase Transition in the Harmonic Model? 

The change from non-Gaussian to Gaussian statistics in our Harmonic model simulations 
is reminiscent of a Kosterlitz-Thouless-Berezinskii transition from a critical region to a 
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Figure 5. Vortices in the Harmonic Model: These figures are snapshots of the 
Harmonic model with N = 3600 at (a) T/J = 1.5, (b) T/J = 2.5, and (c) T/J = 3.5. 
At low temperatures one sees the surprising emergence of a tightly bound vortex/anti- 
vortex pair. The vortex density increases with T and despite a small degree of 
separation of some pairs in (b) and (c), isolated vortices are not observed. The grey 
squares represent a vorticity of -1, the black squares a vorticity of +1. 



paramagnetic phase [44]. The temperature at which Gaussianity is first observed ties in 
with the expression for the susceptibihty derived in [12] 

, = . iV--/*.' (23) 

{a2d ~ 258.6) which has x diverging only for T/J< An. There is, however, no change in 
the behaviour of the magnetization as one passes through this point as one would expect 
for a true phase change and it must be remembered that ()23p is only approximate [12]. 

The KT transition is associated with the onset of topological order, which suggests 
that, if such a transition were to occur in the Harmonic model, a suitable topological 
defect must be found. We have identified one such defect, and, rather surprisingly, it is 
a spin vortex. The definition of a vortex is complicated by the non-periodicity of the 
spins. We assign a 'winding number' to each pair of spins equivalent to the multiple 
of 2iT that must be subtracted from their difference to give a value in the range ±7r. 
The sum over the winding numbers around a closed path is then the vorticity of the 
enclosed region. In a system with = 3600 we observe a tightly bound vortex/anti- 
vortex pair aXT / J = 1.5 (figure Ej). It appears as though the vortices are created by the 
superposition of high energy spin waves, in contradiction to the assumption that spin 
waves and vortices are independent [35,44] (although, as previous studies have generally 
focused on temperatures far below T^t this should have no bearing on their results). 
The vortex density increases with T, eventually becoming so large that it is difficult to 
distinguish individual pairs. At no point do any vortices have a vorticity greater than 1 
and isolated vortices are never observed. 

The lack of vortex pair unbinding can be explained by considering the free energy 
of an isolated vortex. Adapting the standard calculation (see, for example, [45]) to allow 
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for the non-periodic spins, we arrive at an expression for the energy of a vortex, 

E^„, = 2 j7r^(L - a) - JttI In (^^^ , (24) 

where / is a constant of the order of the lattice spacing a. That E'vor diverges as \/iV 
means that it dominates the entropy (which diverges logarithmicaUy with N) at all 
temperatures. Thus isolated vortices are always energetically unstable in the Harmonic 
model. Bound pairs of vortices whose cores are separated by a distance R will have 

^vor-pair ~ 0{R - ]n{R)) (25) 

and are therefore feasible provided R remains small. 

Despite the observation of topological defects, the evidence from statistics and 
the behaviour of the susceptibility, the lack of vortex-pair unbinding must lead to the 
conclusion that no Kosterlitz-Thouless-Berezinskii transition occurs in the Harmonic 
model. 

6. Conclusion 

In this work we have considered the PDF of order parameter fluctuations in the finite 
2dXY model in the spin wave approximation. We have not generally considered the 
effect of spin vortices on the PDF. Bound vortex pairs are irrelevant perturbations in 
the low temperature phase and so might be expected to leave the PDF unaltered, but 
this has not been proved explicitly. For a study of the effect of vortex unbinding on the 
PDF, see [14]. 

The original result of [13] suggested a PDF describing magnetization fluctuations in 
critical finite 2dXY magnets which was truly universal and had been seen to describe the 
statistics of a range of critical systems. We have demonstrated analytically that there 
is in fact a temperature dependence and so it must be concluded that strict universality 
does not hold. However, numerical work shows that the effect of temperature is very 
weak. From a visual point of view the distribution looks very similar across the critical, 
vortex free, region. The temperature dependence established here is weaker than, but 
of the same magnitude as, that claimed by Palma et al, for the full XY model [33]. 
More work is needed to fully understand this difference. 

Despite the weakness of the temperature dependence, we feel that the confirmation 
of its existence is an important result, particularly as much literature has grown up 
around the unusual behaviour of fluctuations in the 2dXY model. Signiflcantly it has 
been shown that the equivalence between the distributions obtained using the cosine and 
linearized order parameters in [13] is the result of the imposition of low temperature in 
both cases - by neglect of MLGs and anharmonic terms respectively. 

We note that the assertion that the Harmonic model is 'vortex free' is true only 
in the sense that vortices do not appear explicitly in the Hamiltonian. However, by 
introducing a deflnition of vortices consistent with inflnitely variable spins, we have 
identifled tightly bound vortex pairs at high temperatures which we conclude are the 
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result of the superposition of high energy spin waves. The crossover to a Gaussian 
statistical regime at T ~ 47r for the Harmonic model, coupled with the fact that the 
susceptibility appears not to diverge above this point, are suggestive of a Kosterlitz- 
Thouless-Berezinskii transition. However, the unbinding of vortex pairs this requires 
is not observed and we argue that isolated vortices in this model are not energetically 
viable. 
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